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Abstract 

The second-order reduced density matrix method (the RDM method) has performed well in 
determining energies and properties of atomic and molecular systems, achieving coupled-cluster 
singles and doubles with perturbative triples (CCSD(T)) accuracy without using the wave-function. 
One question that arises is how well does the RDM method perform with the same conditions that 
result in CCSD(T) accuracy in the strong correlation limit. The simplest and a theoretically 
important model for strongly correlated electronic systems is the Hubbard model. In this paper, 
we establish the utility of the RDM method when employing the P, Q, G, Tl and T2' conditions 
in the two-dimensional Hubbard model case and we conduct a thorough study applying the 4x4 
Hubbard model employing a coefficients. Within the Hubbard Hamiltonian we found that even 
in the intermediate setting, where U/t is between 4 and 10, the P, Q, G, Tl and T2' conditions 
reproduced good ground state energies. 



*Electronic address: maho@riken.jp 
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I. INTRODUCTION 



The second-order reduced density matrix is necessary and sufficient to compute all the 
physical properties that one can compute using the wave-function [1]. Due to its simplicity, 
it has been a dream for quantum chemists to directly determine the second-order reduced 
density matrix instead of using the wave-function, and we believe that it should be simpler 
to determine than solving the Schrodinger equation. 

When an appropriate subset of necessary iV-representability conditions, a term coined by 
Coleman [2] , are used as constraints in a variational calculation of the second-order reduced 
density matrix one is able to compute accurate energies of the second-order reduced matrices 
producing accurate energies and properties. This approach is known as the RDM method 
and has a long history [3, 4]. Unfortunately, the RDM method faded away because no 
algorithm for systematic calculations was available at the time and the iV-representability 
condition was not very well understood. 

After 25 years, in 2001, Nakata et al. formulated the RDM method as the standard form 
of the primal semidefmite programming problem. They performed a systematic study on 
small (few electron) atoms and molecules [5]. They used the P, Q [2], and G conditions [6] 
as the iV-representability constraints that resulted in 120% of correlation energies. These 
promising results led Zhao et al. three years later to include the Tl, and the T2-conditions 
in addition to the P, Q, and G conditions in the RDM method giving results with similar 
accuracy to coupled-cluster singles and doubles with perturbative triples (CCSD(T)) for 
atomic and molecular systems [7-9]. Since then, research along these lines has spread with 
enthusiasm and several papers have been published [10]. 

However, the correlation in molecular systems is not especially strong. We want to 
investigate the robustness of these conditions in predicting accurate energies in the case of 
strong correlation. To test this, we have chosen to employ the Hubbard model [11]. This 
model is interesting not only because of its simplicity, but also its capability of describing 
strong electron correlation. The RDM method has been applied to the Hubbard model 
by Hammond et al. [12], Nakata et al. [9], and Verstichel et al. [13]. Their results very 
accurately described total energies as well as other properties. However, they only treated the 
one-dimensional Hubbard model, which can be solved analytically by the Bethe-Ansatz as 
demonstrated by Lieb and Wu [14] . It can also be treated numerically by the density matrix 
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renormalization group (DMRG) method [16]. As a result, the behavior of the correlation is 
rather well understood [15]. 

The challenge for the condensed matter physics community is, thus, to compute the 
ground state energy and properties of the two-dimensional Hubbard model since no analytic 
results are available as they are in the one-dimensional case. Still, it is an open question, 
but it is believed that two-dimensional Hubbard model is the simplest model that exhibits 
the high-T c superconductivity of copper oxide [17]. The underlying physics of the Hubbard 
Hamiltonian remains a topic of considerable discussion [18]. 

This problem can be reduced to the eigenvalue problem of astronomically large symmetric 
matrices. Extensive numerical studies [19] have been done using the Quantum Monte Carlo 
(QMC) method, the Exact Diagonalization (ED) method (also known as the full configura- 
tion interaction (FCI) method), and the DMRG method. However, we can solve very small 
two-dimensional Hubbard model system without much difficulty To the best of the authors' 
knowledge the largest two-dimensional Hubbard model systems that have been treated are 
the 10 x 10 square lattice by Sorella or 16 x 16 square lattice by Chen et al. [20], 40 to 64 
rectangular or square lattices by the DMRG [21], and the V20 x \^20 by the Exact Diag- 
onalization [22]. Aside from the Exact Diagonalization, the accuracy of the ground state 
energies can be dubious. 

Advantages of the RDM method are: this method calculates the lower bound to the 
FCI energy in the same basis set whereas all of the other methods give upper bounds, 
thus this approach is complementary to the former methods. This method does not require 
extrapolation to the absolute zero-temperature. It does not suffer from minus sign problem 
in QMC [30]. It does not depend on the choice of lattice which may appear in DMRG 
calculation [31]. 

In this paper, we calculated the total energies of the two-dimensional Hubbard model 
using the RDM method and compared them to the exact results from ALPS [24] to examine 
whether the P, Q, G, Tl and T2' conditions are physically important in strongly correlated 
system. The rest of paper is organized as follows. In Section II, we briefly review the RDM 
method, the iV-representability conditions, semidefinite programming, and the Hubbard 
models. The results and discussion are shown in Section III. The conclusions are in Section 
IV. 
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II. THEORY 



A. Reduced density matrices 

The second-order density matrix is an example of a broader class of density matrices. 
The most general form is the M-th order density matrix. This has the form 

(M) T hi2-i„ J_ATr| a \ a \ . . .J „. . . . • fl . I \IA 

1 nn-m ^|\^l u n u « 2 u i m u 3m u n a n\^ I 

The second-order reduced density matrix is an important special case. This is the reduced 
density matrix that we are utilizing. Explicitly, it has the form: 

When dealing with the 1-body terms (present in most Hamiltonians of interest and several 
properties operators), the second-order reduced density matrix reduces to the first-order 
reduced density matrix defined as: 

7} = <*|al % |*>, 

where a* and a denote the creation and annihilation operators, respectively, and \& is the 
TV-particle antisymmetric wave-function. Note that it is usually denoted by 7* instead of by 

rj. 

The second-order reduced density matrix has seen renewed interest for computing dy- 
namical properties of a quantum mechanical system governed by the electronic Hamilto- 
nian. When this descriptor was first introduced as a descriptor for electronic structure it 
was met with enthusiasm [3] . Unfortunately, when the RDM method was applied to nuclear 
systems like 24 Mg, 28 Si, the energies were found to be far below the expected value [4]. This 
is because the second-order reduced density matrix that resulted from these calculations 
did not originate from any wave-function [2] ! Every reduced density matrix of interest must 
result from some wave-function (this wave- function is known as the ancestor wave- function). 
The problem of reduced density matrices not arising from ancestor wave-functions is what 
Coleman [2] coined the iV-representability problem. Currently the necessary and sufficient 
conditions that guarantee iV-representability are not known in any practical form [23] . For- 
tunately, several necessary conditions are known. Using only the P, Q, G, Tl, and T2' 
(necessary) conditions have been shown to reliably obtain chemical accuracy [7, 9]. 



B. iV-representability conditions 



7V-representability is the necessary and sufficient conditions that a density matrix orig- 
inates from some (ancestor) wave-function [2]. For the first-order density matrix to be N- 
representable its eigenvalues should lie in the closed interval [0, 1] [2, 25]. Since we know these 
conditions in an implementable form and due to Gilbert's theorem [26], one can construct 
a method using only the 1-RDM. This method is sometimes referred as the density-matrix 
functional theory (DMFT) method [33]. The iV-representability conditions are not limited 
to density matrices. The iV-representability conditions for the electron density are known 
[26]. For the wave- function itself they are very simple, simply ensure the basis functions are 
square integrable and antisymmetric (change sign) with respect to the interchange of any 
two electron coordinates (Pauli principle). 

Unfortunately, the second-order density matrix iV-representability conditions are not 
known in any useful form (i.e. an uncountable set of conditions) [6, 23]. However, many 
necessary conditions are known. Some trivial conditions are trace conditions; 

E^ = Ar ' £rg = Jv(Jv-i)/2, (2.i) 

i ij 

and 

k 

An incomplete list of necessary conditions alone are not enough to guarantee N- 
representability. However, using them within the RDM method gives strict lower bounds to 
the energy. The general strategy within the RDM method is to choose necessary conditions 
that are easily implementable, computationally inexpensive, and result in accurate energies. 
Of course, the more necessary conditions used the better the answer (though how much the 
energy is improved depends on the system being investigated and the necessary condition 
being used). The most commonly utilized conditions used within the RDM method are 
positive-semidefinite type of iV-representability conditions; the P, Q [2] and G conditions 
[6]. 

The P-condition is formulated by starting from the simple fact that if A is an arbitrary 
one-particle operator, then the expectation value of A* A should be non-negative, 

(A* A) = Tr(A^A)T > 0. 
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If we restrict A to A = Cijdidj, for an arbitrary set of real numbers c^, then 

(AU> = X) CyC „rg>o. 

should be satisfied. Therefore, T should be positive semidefinite. Explicitly, the P-condition 
is: 

where we used X >z O to indicate it is positive semidefinite. If we restrict A to A = 
^Cijoloj, then likewise, the Q-condition [2] is explicitly: 

The Q matrix should also be positive semidefinite. The G-condition can be derived if we 
restrict A to A = Y2 c ij a l a j- The Tl- and T2-conditions are derived in a somewhat more 
involved way [7]. If we take A as A = J2ijk c ijkCi m n(a\ajala n a m a£ + a n a m ae,a\ a ) a k)i which is 
a positive semidefinite three-particle operator, but it cancels out the genuine three-particle 
part. Thus we can evaluate it using the second-order reduced density matrix. This is the 
Tl condition. In the same way, if we take A = Y^ijk c ^ c irnn{a\a\a k a\a m at + d n a m aAAj a k) , 
which is also a three-particle operator, and then cancels the three-particle part again, then 
we have the T2-condition. Usually, the T2' condition is used instead of the T2-condition. 
The T2'-condition is an enhancement of the T2-condition that arises from the addition of 
the one-particle operator [27]. 

C. The RDM Method 

As stated in the previous section, the necessary and sufficient conditions for N- 
representability are not known in any useful form for the second-order reduced density 
matrix. However, many necessary conditions are known. Selecting necessary conditions for 
7V-representability would be how we find the set S^, the set of approximately (necessary) 
7V-representable second-order reduced density matrices. Then the RDM method for the 
ground state is the minimization of the total energy subject to: 

Eg = min {Tri/r} , 

= min K>$7j+ £ w^T^l, 

N I ij i\iihh ) 
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where we used the definition of the (electronic) Hamiltonian as follows: 

ij i\iihh 

Finding and implementing conditions that expedite calculations is the strategy for increasing 
the accuracy of the RDM method. One can formulate a method for each of the M-th order 
reduced density matrices if the appropriate boundary conditions are included. All of the 
matrices used in the RDM method are positive semidefinite. As a result of this, a semidefinite 
optimization program is used [5, 7]. 

As mentioned in the previous section, an interesting difference between the variational 
methods and the RDM method are we always obtain lower bounds to the exact energy 
whereas the variational methods gives the upper bounds [1, 2]. 

D. Semidefinite Programming 

Semidefinite programming (SDP) is a convex optimization problem, a smooth generaliza- 
tion of the linear programming. The linear programming optimizes a linear functional of a 
non-negative vector whereas the SDP optimizes a linear functional of a non-negative, (i.e., 
positive semidefinite symmetric) matrix. An SDP problem can be defined as: 

(P) min : Y% =l c k x k 

s.t. : X = T^ =1 F k x k -F ,Xy O, 

where the symbol S n is the space ofnxn symmetric matrices. We use X >z 0(X >- O) 
to indicate X e S n is positive semidefinite (positive definite, respectively). The Lagrangian 
dual (D) of the problem (P) can be derived as: 

(D) max : F • Y 

s.t. : F k • Y = c k (k = 1, . . . , m) 
YhO. 

Here we used the Hilbert-Schmidt inner product U • V defined as Yli=i Y^=i UijVij for U 
and V in § n . 

To solve the (P) formulation or the (D) formulation, we employ the primal-dual path- 
following interior-point method (PDIPM). This algorithm solves these two problems, (P) 
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and (D), simultaneously in polynomial time. It is a widely accepted method because there 
are also many efficient implementations [28]. One of the most efficient implementation 
is SDPARA developed Fujisawa et al. [29]. We used this program to solve large scale 
semidefinite programming problems arising from condensed matter physics. 
Here we briefly sketch a framework of the PDIPM: 

Step 0: Choose an initial point x°, X°, Y° with X° y O,Y y O. Set h = and choose 
the parameter 7 e (0, 1). 

Step 1: Evaluate the Shur Complement Matrix B G S n by the formula 

B, l] = {{X h )- l F l Y h ).F r (2.3) 

Step 2: Solve the linear equation Bdx = r. Using its solution dx, compute dX, dY and 
obtain the search direction (dx,dX,dY). 

Step 3: Compute the maximum step length a to keep the positive semidefiniteness; a = 
max{a G [0, 1] : X h + adX y 0,Y h + adY y O}. 

Step 4: Update the current point by (x h+1 , X h+1 , Y h+1 ) = (x h , X h , Y h )+^a(dx, dX, dY). 

Step 5: If (x h+1 , X h+1 ,Y h+1 ) satisfies the stopping criteria, output it as a solution. Oth- 
erwise, set h — h + 1 and return to Step 1. 

Since basically, all the iV-representability conditions can be written as linear inequalities 
[6], the RDM method with any kind of iV-representability condition can be formulated as 
semidefinite programming problem. 

There are two ways of formulating the RDM method as a standard type semidefinite pro- 
gramming problem. One is the primal formulation [5] and the other is the dual formulation 
[7] . The primal formulation is somewhat more involved than the dual formulation, however, 
the number of variables is reduced considerably. In either case, the formulation would result 
in large size SDPs. 

E. The Hubbard Model 

The Hubbard model is a simple lattice model that was formulated to model strong correla- 
tion [11]. It was used initially to show the behavior of d electrons. In its original application 
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the Hubbard model was used to describe electrons in solids. It is now used for predicting 
superconductivity, particularly, it is used to describe the transition between conducting and 
insulating bands [15]. The Hubbard Hamilton for fermions has the form: 



where U and t are real parameters, L is the number of sites on the lattice, and means 
summing over every i-th and j-th sites that is the nearest neighbors. The Hamiltonian is 
made up of two terms. The first term is the kinetic term allowing for tunneling. It is also 
referred to as the hopping term. The second term is the on-sight repulsion term. The tight 
binding model is a special case of the Hubbard model. It arises when the on-site repulsion 
term is neglected (U — 0). It has this name because it describes tightly bonded electrons in 
solids. There are limited interactions of each particle with its neighbors. This is similar to 
orbitals within a free atom in the standard linear combination of atomic orbitals (LCAO) 
model. Alternatively, if t — we obtain the one-site model. As the name implies, all the 
sites are independent of each other (the neighbors no longer matter). Each site may contain 
electrons, 1 alpha electron, 1 beta electron, or 1 alpha and 1 beta electron. When t ^ 
then only the ratio of U/t matters [15]. It is a model that is easily adaptable to using the 
RDM method. 

III. RESULTS AND DISCUSSION 

We targeted the repulsive 2-dimensional 4x4 Hubbard model with periodic boundary 
condition, and using hopping term that is non-zero for nearest neighbors and zero otherwise. 
The number of electrons in the system is 16 (half-filled) and the total spin of the system 
being S 2 = 0. Then, we solved for energies using Hamiltonians with U/t values ranging from 
0.01 to 100. Note that in the bipartite lattice case, there is a symmetry in U/t; changing 
the sign of t does not alter the physics (in that only the ratio of U/t matters), thus the same 
2-RDM is obtained since the gauge transformation on the one sub lattice A is: 



can change t to —t. Thus, in this paper, we fixed to t — 1. We employed the P, Q, G, Tl, 
and T2' as A^-representability conditions and three types of combinations with trivial N- 



L 



L 




Wi,a}A -> {-<,*} A 
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representability conditions (eqs. 2.1 and 2.2) are examined; (i) the P, Q and G conditions, 
(ii) the P, Q, G, and Tl , (iii) the P, Q, G, Tl, and T2' conditions. To compare the accuracy 
of each calculation we also performed the Exact Diagonalization method using ALPS [24]. 

In Table I we show the details of the cluster machine used to solve the large SDPs of the 
2-dimensional Hubbard model, we show the size of the problem in Table II, and we show 
typical time elapsed to solve the problem in Table III. As we see, the number of constraints 
do not change in the dual formulation [7], however the number of variable of the matrix 
becomes larger when we add the T2' conditions; the order of size is the same as in the third- 
order matrices. No simple collapse can be found as the Tl condition's case. The elapsed 
time shows that to solve using the P, Q, and G-conditions that it took approximately 1000 
to 2000 seconds, to solve using the P, Q, G, and Tl-conditions it took approximately 4000 to 
5000 seconds, and to solve using the P, Q, G, Tl and T2'-conditions it took approximately 
20, 000 to 30, 000 seconds. Since the convergence criteria is well understood in the PDIPM, 
the number of iterations required to converge is dependent of the size of the problem, and 
not on the form of the Hamiltonian. This is another feature of the RDM method when 
solved using the PDIPM. 

We have omitted the results here, but the numerical quality of the solutions (the con- 
vergence) is approximately eight significant decimal digits. Theoretically, the primal-dual 
gap, primal and dual feasibility of the SDP (P) and (D) should be zero or numerically they 
should be sufficiently small (as determined by convergence criteria) [29] . These values for our 
solutions are very small in our results; typical values for the primal-dual gap are: 1.0 x 10~ 10 
to 1.0 x 10 -15 , and for the primal and dual feasibility are: 1.0 x 10~ 7 to 1.0 x 10~ 9 . 

In Table IV, we show the total energy of various values of U/t by (i) the ED method, (ii) 
by the P, Q, and G, (iii) by the P, Q, G, and Tl, and (iv) by the P, Q, G, Tl, and T2' 
conditions. The difference in the energies from the ED method, AEpqg = Epqg — Eed, 
AE PQGT1 = EpQ GTX - E ED , and AE PQGT1T2 > = E PQGT1 - E ED are shown. Since adding 
7V-representability conditions would result in increasing (and thus improving) the energy, 
therefore, Epqg < Epqgti < Epqgtitv < Eed always holds theoretically, and we confirm 
this result numerically. 

In the weak coupling limit \U/t\ « 1, same as in the one-dimensional case, the total 
energies calculated by the RDM method almost coincide with the ED method; at U/t = 0.01, 
the total energies from the RDM method is almost the same as the ED method. As the 
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TABLE I: The following machines and SDP solver have been used in calculating the ground state 
energy of the Hubbard model. 



Cores 


16 Nodes, 32 CPUs, 128 CPU cores 


CPU 


Intel Xeon 5460 3.16GHz (quad cores) x 2 / node 


Memory 


48GB / node 


NIC 


GbE x 2 and Myrinet-lOG x 1 / node 


OS 


CentOS 5.8 for x86_64 


SDP Solver 


SDPARA 7.3.2 RC2 



coupling becomes larger, the total energy decreases rather quickly for E P q G , especially when 
U/t is larger than I. The worst energies are obtained when 4 < U/t < 10. By adding the 
Tl condition, the absolute difference in energy reduces by a factor of two for most values of 
U/t. The most accurate results were obtained by when using the P, Q, G, Tl and T2' as 
the iV-representability conditions. The largest error was obtained at U/t — 8 (—2.39 x KU 1 ) 
or, equivalently, —1.8 x 10~ 2 per site. This is still a very good energy [19]. In the range 
where the worst energies were obtained, 4 < U/t < 10, and is an important range because 
in this parameter region quantum phase transition from metal (tight-binding one) to Mott- 
insulator phase (and to the Neel phase) typically occurs. For this reason, the calculation of 
intermediate coupling is difficult [15]. 

We cannot show the results for U/t > 1000, since we faced the numerical difficulty in 
this region. In the high correlation limit (\U/t\ — > oo, all the states are nearly degenerated. 
In this case, SDP problems also become degenerated. This results numerical instability and 
we need to perform high precision calculations [9]. 

It is not clear what happens if we add more sites and electrons. The RDM method with 
P, Q, and G is not size-extensive; the total energy should scale if the size of the system is 
scaled [32]. In molecular cases, reasonable bounds were found, but there is no guarantee 
to have a lower bound on the energy of the unit lattice. Size-extensively may be recovered 
when we add unitary invariant iV-representability condition [34]. 
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TABLE II: Size of the SDPs for the 4x4 Hubbard model using various iV-representability condi- 
tions. The number of constraints, the size of variable matrix, and the number of block matrices in 
the variable matrix and the maximum size of the block in the variable matrix are shown. 



iV-representability 


constraints 


size of variables matrix 


# of blocks 


maximum block 


P, Q, G 


47688 


2634 


14 


512 


P, Q, G, Tl 


47688 


7594 


18 


1920 


P, Q, G, Tl, T2' 


47688 


23498 


22 


6032 



TABLE III: Typical elapsed times to solve the SDP problems. 



iV-representability 


Time (s) 


P, Q, G 


1000 ~ 2000 


P, Q, G, Tl 


4000 ~ 5000 


P, Q, G, Tl, T2' 


20, 000 ~ 30, 000 



IV. CONCLUSION 

The P, Q, G, Tl and T2' conditions were found to give reliable accuracy for the strongly 
correlated electron case as espoused by the two-dimensional Hubbard model. The most 
difficult case for these conditions to resolve was the intermediate case where 4 < U/t < 10. 
The accuracy was more than satisfactory when using the P, Q, and G conditions. The 
deviation of the energy is 0.1 unit for one lattice in this region. By adding the Tl condition, 
the deviations from the exact energies were improved by approximately a factor two from 
using only the P, Q and G-conditions. Prominent results were obtained further adding 
T2' condition. The deviations from the exact energies are approximately 0.01 unit for 
intermediate couplings. If one could formulate an A^-representability condition for a tightly 
bound electron with some small probability of mobility this would increase the accuracy 
further. Nevertheless, the results in these cases are satisfactory and better in all other 
cases. Since solving semidefinite programming problems from condensed matter physics 
becomes extremely large, then reducing the size by employing symmetry of the system 
is very important. The size of systems we calculated are too small to be conclusive and 
still not yet competitive with other methods. Anyway, it is a very good challenge for 
the optimization community and high performance computing community. The reduced- 
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TABLE IV: The total energy of the ground state of the 2-dimensional Hubbard model with various 
U/t and various methods: Eed is by the ED method, Epqg is by the RDM method using the P, Q, 
and G-conditions, , Epqgti is by the RDM method using the P, Q, G, Tl-conditions, , Epqgtit2' 
is by the RDM method using the P, Q, G, Tl and T2'-conditions. AEpqg = Eed - Epqg, and 
so on, and these AEpqg, AEpqgti and AEpqgtit2' should be negative. 



U/t 




EpQG 


EpQQXl 


E POHTA TO' 


AEp QG 


AEpQGTl 


AE P QGT1T2' 


0.01 


-23.9656 


-23.9657 


-23.9657 


-23.9657 


-2.39 x 10" 


•5 


-1.59 x 10" 


•5 


-1 X 10~ 7 




0.1 


-23.6587 


-23.6606 


-23.6599 


-23.6587 


-1.98 x 10" 


■3 


-1.22 x 10" 


•3 


-1 X 10~ 5 




0.2 


-23.3221 


-23.3298 


-23.3268 


-23.3221 


-7.74 x 10" 


•3 


-4.74 x 10" 


■3 


-1 x 10~ 5 




0.5 


-22.3402 


-22.3858 


-22.3682 


-22.3411 


-4.55 x 10" 


■2 


-2.79 x 10" 


•2 


-8.64 x 10" 


■4 


0.8 


-21.3991 


-21.5090 


-21.4666 


-21.4024 


-1.10 x 10" 


■1 


-6.75 x 10" 


•2 


-3.31 x 10" 


■3 


1 


-20.7936 


-20.9584 


-20.8953 


-20.7998 


-1.65 x 10" 


•1 


-1.02 x 10" 


•1 


-6.13 x 10" 


•3 


2 


-18.0176 


-18.5478 


-18.3522 


-18.0535 


-5.30 x 10" 


•1 


-3.35 x 10" 


•1 


-3.60 x 10- 


■2 


3 


-15.6367 


-16.5790 


-16.2473 


-15.7243 


-9.42 x 10" 


•1 


-6.11 x 10" 


•1 


-8.77 x 10" 


•2 


4 


-13.6219 


-14.9454 


-14.4941 


-13.7711 


-1.32 




-8.72 x 10" 


•1 


-1.49 x 10" 


■1 


5 


-11.9405 


-13.5745 


-13.0214 


-12.1479 


-1.63 




-1.08 




-2.07 x 10" 


■1 


6 


-10.5522 


-12.4134 


-11.7728 


-10.8045 


-1.86 




-1.22 




-2.52 x 10" 


•1 


7 


-9.41048 


-11.4208 


-10.7033 


-9.69084 


-2.01 




-1.29 




-2.80 x 10" 


■1 


8 


-8.46888 


-10.5641 


-9.77809 


-8.76192 


-2.10 




-1.31 




-2.93 x 10" 


■1 


9 


-7.68624 


-9.81887 


-8.97982 


-7.98034 


-2.13 




-1.29 




-2.94 x 10" 


•1 


10 


-7.02900 


-9.16556 


-8.28788 


-7.31630 


-2.14 




-1.26 




-2.87 x 10" 


■1 


100 


-7.68192 x 10" 1 


-1.21923 


-8.75302 x KT 1 


-7.83706 x 10- 1 


-4.51 x 10" 


•1 


-1.07 x 10" 


•1 


-1.55 x 10" 


■2 



density matrix method has been a promising method of quantum chemistry, however, in 
this paper, we showed that this method is also promising for condensed matter physics 
where the electron correlation are very strong. 
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